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SADDLEPOINT  APPROXIMATIONS  FOR  VARIOUS  STATISTICS  OF  DEPENDENT 
NON— GAUSSIAN  RANDOM  VARIABLES:  APPLICATIONS  TO  THE  MAXIMUM 

VARIATE  AND  THE  RANGE  VARIATE 

INTRODUCTION 

Evaluation  of  the  joint  cumulative  distribution  function 
(CDF)  or  joint  exceedance  distribution  function  ( EDF )  from  the 
corresponding  joint  probability  density  function  (PDF)  is  a  very 
difficult  task  in  M  dimensions  when  M  is  approximately  greater 
than  four.  Even  in  the  simplest  nontrivial  case,  namely,  when 
the  joint  PDF  is  Gaussian  with  known  non-diagonal  covariance 
matrix  and  mean  vector,  the  required  integrals  over  an  M- 
dimensional  subspace  are  not  amenable  to  numerical  integration. 

However,  a  number  of  probabilistic  problems  are  amenable  to 
the  analytic  evaluation  of  the  joint  moment  generating  function 
(MGF)  of  the  M  random  variables  (RVs)  of  interest,  even  though 
the  corresponding  joint  PDF  is  not  available  in  any  useful 
analytic  form  (reference  1).  In  these  cases,  it  would  be  very 
useful  to  have  a  means  for  determining  the  joint  CDF,  joint  EDF, 
and/or  joint  moments  directly  from  the  M-dimensional  joint  MGF, 
without  having  to  bother  with  the  intermediate  joint  PDF.  This 
report  presents  just  such  a  technique  for  accomplishing  this 
goal  by  extending  Parseval's  theorem,  relating  one-dimensional 
Fourier  transforms,  to  M-dimensional  Laplace  transforms  with 
movable  Bromwich  contours. 
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The  additional  freedom  afforded  by  moving  contours  in  complex 
M-dimensional  space  allows  for  use  of  an  M-dimensional 
saddlepoint  (SP)  of  the  integrand  relating  the  joint  MGF  to  the 
statistic(s)  of  interest.  A  saddlepoint  approximation  (SPA)  can 
then  be  developed  about  this  special  point  for  the  desired 
statistical  quantity  of  interest.  Although  this  procedure 
requires  the  numerical  determination  of  an  M-dimensional  SP,  it 
affords  a  reasonable  practical  means  for  approximately  evaluating 
some  joint  probabilities  and  joint  moments  that  have  been 
previously  inaccessible.  An  excellent  review  of  the  SP  method  is 
given  in  reference  2,  in  addition  to  a  thorough  bibliography. 

The  calculation  of  these  M-dimensional  statistics  is  not 
bought  cheaply.  It  is  sometimes  required  that  multiple 
M-dimensional  SPs  be  determined  to  evaluate  an  EDF  or  PDF  at  a 
single  point  in  probability  space.  Also,  the  storage 
requirements  and/or  execution  times  can  rapidly  grow  unmanageable 
as  the  number  of  dimensions,  M,  increases.  Finally,  care  must  be 
taken  to  ensure  that  the  appropriate  SP  is  located  in  the  correct 
region  of  M-dimensional  space;  this  allowed  region  of  analyticity 
varies  with  the  particular  statistic  under  investigation. 
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PROBABILITY  NOTATION 


Random  variables  will  be  denoted  by  boldface  type.  Thus, 
real  random  vector  (RV)  z  is  an  Mxl  column  vector 

z  =  [z1  • • •  zM] '  ,  (1) 

with  random  components  {zm},  m=l:M,  that  may  be  statistically 
dependent  on  each  other.  On  the  other  hand,  z  =  [z1  •••  zM ]'  is 
an  Mxl  vector  of  ordinary  real  variables  {z^}.  The  joint  PDF  of 
RV  z  at  argument  z  is  denoted  as 

Pz(z)  “  Pz( zi • • ♦ • • ZM)  •  (2) 


The  corresponding  joint  CDF  of  RV  z  at  argument  z  is 


JM 


cz(z)  -  cz(Sl . ZM)  -  J  dUj  ...  J  duM  Pz(Uj . UM) 


-  <  . . H  <  ZM>  ' 


(3) 


where  Pr(e)  is  the  probability  of  event  e.  The  corresponding 
joint  EDF  of  RV  z  at  argument  z  is 


ez(Z>  -  ez(zl . ZM>  -  J  dul  •••  J  dUM  Pz(Ul . V  * 


M 


"  Pr(z!  >  >  ZM}  * 


(4) 


It  should  be  noted  that  c  (z)+e  (z)  <1,  generally,  in  M 

z  z  * 

dimensions  for  M  >  1. 
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MIXED  PROBABILITY  FUNCTIONS 

It  will  be  necessary  later  to  introduce  and  utilize  some 
mixed  probability  functions.  For  example,  the  following  quantity 
is  M-l  parts  CDF  and  1  part  EDF : 

Pr ( 21  <  zl'*“'  ZM-1  <  ZM-1 '  ZM  >  ZM)  ” 

Z1  ZM-1 

*  J  dUl  •"  J  duM-l  J  duM  Pz(ul . V  •  <5> 

—00  —  00  9 

M 

On  the  other  hand,  the  following  quantity  is  1  part  PDF  and  M-l 
parts  CDF: 

z2  ZM 

J  du2  •••  J  duM  Pz(zru2 . UM)  .  (6) 

—  00  —oo 


Finally,  although  the  following  quantity  is  not  encountered  in 
this  report,  it  is  listed  to  demonstrate  the  great  generality  of 
the  probability  procedure  to  be  presented  here,  namely. 


z^  ® 


J  dul  J  du4  J  du5  J  du7  J  du9  PZ(U1'Z2'Z3'U4'U5'Z6'U7'Z8'U9'Z10) 

(7) 


2  4  z5 


—  00  —00 


is  five  parts  PDF,  three  parts  CDF,  and  two  parts  EDF.  It  should 
be  observed  that  all  of  the  probabilistic  quantities' here  are 
special  cases  of  the  M-dimensional  integral 


J 


du. 


I 


du 


M 


PZ(UX' 


'V 


9(u, 


•'V 


-J 


du 


pz(u) 


g(u).  (8) 
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TRANSFORM  DOMAIN 


The  transform  domain  is  denoted  in  terms  of  Mxl  complex 
vector  X  =  [X1  •••  XMJ'.  In  particular,  the  joint  MGF  ^Z(X) 
corresponding  to  joint  PDF  pz(z)  is  given  by  expectation 


M  N 

S.  x»> 


fi  (X)  =  E  exp(X'  z)  =  E  exp| 
z 

=  J  dzx  •••  J  dzM  pz(z1, ...,zM)  exp( X'  z) 


(9) 


When  the  imaginary  part,  X^,  of  vector  X  is  zero,  integral  (9) 

exists  for  Re(X)  e  R  ,  where  M-dimensional  real  region  R  always 

//  " 

includes  the  origin  X  <=  0.  It  immediately  follows  from  the  form 
of  equation  (9)  that  joint  MGF  //z(X)  will  exist  for  all  values  of 
vector  \i  when  Re(X)  s  R^.  This  M-dimensional  X  region  of 
definition  of  joint  MGF  u  (X)  is  called  the  region  of 
analyticity,  ROA that  is,  ROA(/uz)  is  the  M-dimensional  set 
of  complex  X  values  such  that  Re(X)  e  R^. 


The  joint  PDF  p  (z)  at  vector  argument  z  may  be  obtained  from 
z 

the  joint  MGF  //  (X)  by  means  of  M-dimensional  inverse  Laplace 
z 

transform 


PZU)  = 


( i2n) 


M 


J 


dX, 


I 


dXM  /yz(Xi,  •  •  •  ,XM)  exp(-X'z)  ,  (10) 


'M 


where  Bromwich  contours  {Cm},  m=l:M,  initially  lie  in  the 
ROA ( /w 2 )  •  The  m-th  Bromwich  contour  Cm  parallels  the  imaginary 
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axis  Xmi  and  remains  in  the  ROA(yuz).  However,  the  freedom  to 
move  all  the  contours  {Cm)  within  the  ROA  allows  for  an 
advantageous  choice  of  locations,  namely,  through  the 
M-dimensional  real  SP  of  the  integrand  of  equation  (10)  that  lies 
in  ROA(//z);  see  reference  1  for  an  example. 

The  corresponding  joint  cumulant  generating  function  (CGF)  of 
RV  z  is  given  by 

XZ(X)  =  In  /t/z(X)  for  X  e  ROA(/;z)  .  (11) 
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REAL  NONLINEAR  TRANSFORMATIONS 


Suppose  random  vector  z  is  subjected  to  the  real . nonlinear 
transformation  g(u)  =  g(u^,...,uM)  according  to 

y  =  g(sir...,*M)  »  (12) 

leading  to  real  scalar  RV  y.  Two  examples  are 

y  =  j-j  +  •••  +  z^j  and  y  =  max(  z^ , . . .  /  zjj)  •  (13) 

The  average  value  of  RV  y  in  equation  (12)  is  given  by 
a  =  E(y)  =  E  g( z^ , . . . , z^)  = 

“  J  dui  •**  J  duM  pz(u1,...,uM)  g(u1,...,uM)  .  (14) 

This  is  the  same  M-dimensional  integral  encountered  in  equations 
(5)  through  (8).  However,  M-dimensional  integral  (14)  is 
virtually  always  too  difficult  to  determine  analytically;  also, 
for  large  M,  it  is  extremely  difficult  to  evaluate  numerically 
due  to  storage  and  execution  time  limitations.  Accordingly,  an 
alternative,  more  useful  form  for  integral  (14)  will  be  developed 

Substitute  expression  (10)  for  joint  PDF  pz  into  equation 
(14)  and  interchange  the  M— dimensional  integrals  to  obtain 

*  -  ,  .  :".M  J  dXl  •”  J  dXM  "z<Xl . V  X 

<i2n)  J  c 

1  M 

x  J  dux  •••  J  duM  g(u1,...,uM)  exp(-X'u)  .  (15) 
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The  outer  integrals  in  equation  (15)  require  that  X  be  kept  in 
the  ROA  of  // z ( X )  -  At  the  same  time,  the  inner  integrals  on  u  in 
equation  (15),  to  be  denoted  by  gamma  function  y(X),  will 
converge  only  for  vector  X  within  a  restricted  region  of 
M-dimensional  space.  That  is,  gamma  function 

y(X1, . . . ,XM)  =  J  duj  •••  J  duM  g(u1,...,uM)  exp(-X'u)  (16) 

exists  only  for  X  e  ROA(y).  Use  of  relation  (16)  in  equation 
(15)  yields  the  alternative  expression  of  interest  for  average  a: 


‘  (i2n)«  l 


dX. 


I 


dX^j  // g  (  X^  ,  . . . ,  Xj^ )  y  (  X^  ,  ♦  •  • ,  Xjj  ) 


(17) 


'M 


In  more  compact  notation,  from  equations  (14)  and  (17), 

a  *  f  du  Pz(u)  9(u)  =  — - — jr  f  dX  //_(  X)  y(  X)  ,  (18) 

J  (i2n)W  J  2 

where  M-dimensional  contour  C  must  lie  in  the  intersection  of  the 

ROAs  of  //  (X)  and  y(X);  that  is,  C  e  ROA ( //  , y )  s  roa(//  )  n  ROA(y). 
"  z  z 


The  relation  (18)  is  exact;  there  are  no  approximations 

involved  in  this  identity.  However,  the  M-dimensional  integral 

on  X  will  generally  not  be  capable  of  evaluation  analytically 

either;  but  the  fact  that  the  M-dimensional  contour  C  in  equation 

(18)  can  be  moved  about  in  the  restricted  ROA (//  ,y)  affords  the 

z 

chance  of  developing  an  SPA  for  the  last  term  in  equation  (18) 
about  a  SP  in  this  region. 


8 


PROS  AND  CONS  OF  EQUATION  (18) 

The  first  integral  in  equation  (18),  1^,  has  an  immediate 
physical  interpretation  as  the  average  of  nonlinear  transforma¬ 
tion  g,  whereas  the  second  integral  in  equation  (18),  has  no 

direct  physical  interpretation.  Also,  it  is  easy  to  directly 
specify  nonlinearity  g  in  1^,  whereas  gamma  function  y  in  I2  must 
be  obtained  via  the  M-dimensional  Laplace  transformation  (16). 

On  the  other  hand,  many  probabilistic  problems  have  no 

convenient  expression  for  the  joint  PDF  p  ,  whereas  the  joint  MGF 

z 

fj  can  sometimes  be  obtained  in  closed  form;  see  reference  1,  for 
z 

example,  where  the  joint  MGF  for  M  quadratic  and  linear  forms  in 
K  dependent  Gaussian  RVs  is  derived.  Also,  whereas  development 
of  an  SPA  for  in  equation  (18)  is  difficult  due  to  the 
discontinuities  inherent  in  typical  g  functions,  the  SPA  for  I2 
in  equation  (18)  can  be  developed  more  easily  by  moving 
M-dimensional  contour  C  to  an  appropriate  SP  in  X  space. 

The  availability  of  a  closed  form  for  joint  MGF  jj  , 

Z 

combined  with  the  lack  of  knowledge  of  joint  PDF  p  ,  weighs  very 
heavily  for  form  I2  over  1^  in  equation  (18).  The  major  drawback 
with  I2  is  in  getting  the  M-dimensional  Laplace  transform  y  of 
nonlinearity  g;  however,  for  an  important  class  of  physically 
meaningful  problems,  the  g  function  is  separable  in  its  arguments 
u  =  [u^  •••  uMlf,  thereby  allowing  reduction  of  equation  (16)  to 
M  one-dimensional  Laplace  transforms,  a  very  workable  approach. 
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EXAMPLES  OF  TRANSFORMATION  g(u) 


Example  1.  9Z (^ , . . .  ,uM)  =  5(u1-z1)  •••  $(uM~zM)  ,  (19) 

where  vector  z  =  ^  •••  z MJ '  represents  an  arbitrary  point  in 
real  M-dimensional  space.  Substitution  of  equation  (19)  in  the 
left-hand  side  of  equation  (18)  yields 

a  =  J  du  pz(u)  gz(u)  =  P2(z1,...,zM)  ,  (20) 


which  is  simply  the  joint  PDF  of  RV  z  at  argument  z.  At  the  same 
time,  substitution  of  equation  (19)  into  equation  (16)  yields 


r  (X)  =  exp(-X'z)  =  exp 


M 


-  n 

.  m=l 


X  z 
m  m 


for  any  { Xm > 


(21) 


That  is,  the  ROA(y)  is  infinite  for  this  particular  impulsive  g 
example  in  equation  (19);  therefore  R0A(a/z,y)  =  ROA (//2).  Use  ol 
equation  (21)  in  equation  (18)  yields  the  alternative 


a 


M  I  dX  ^z(X)  exP(-Xf2)  / 
(i2n)  l  2 


(22) 


which  is  the  usual  expression  (10)  for  the  joint  PDF  p  (z). 


Nothing  new  has  been  learned  from  this  particular  example  for 
nonlinearity  g.  However,  the  method  developed  here  will  be  used 
for  all  the  following  examples,  without  additional  explanation. 

A  key  point  is  to  note  the  ROA(y)  when  equation  (16)  is  evaluated 
for  the  particular  g  example  under  consideration.  This  must  be 
coupled  with  ROA(//z)  according  to  ROA(//z,Y)  =  ROA Uz)  n  ROA(y). 
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Example  2.  gz(u1, . . . ,uM)  =  U(z1-u1)  U(zM~uM)  , 


where  unit  step  function  U(x)  is  1  for  x  >  0  and  is  0  for  x  <  0. 
Then,  average  (14)  becomes 


■  J  dUl  "•  J  dUM  PZ(U1 . V  -  CZ(Z1 . H>  ' 


which  is  the  joint  CDF  of  RV  z  at  argument  z  =  [z^  •••  zMJ'.  The 
gamma  function  (16),  corresponding  to  nonlinearity  (23),  is 


r2(X)  -  J  dUj  •••  J  duB  exp(-X1  Ul  -  •••  -XM  uH)  = 


exp(-X1  z1) 


exp(-X  z  )  /X-.  .  . 

- - exp(-X'z)/  rr<-xm) 


if  Re(X  )  <  0  for  m=l:M.  Therefore,  average 


a  =  c 


,(z)  =  — - — j7  [  dX  u  (X)  exp(-Xfz)/  f-f(-X  )  , 
( i2n )M  l  /  m=l  m 


provided  that  the  m-th  Bromwich  contour  C  stays  to  the  left  of 
the  origin  for  m=l:M,  and  also  stays  within  the  ROA  of  joint  MGF 
//Z(X).  Relation  (26)  provides  a  method  for  obtaining  the  joint 
CDF  cg(z)  at  arbitrary  M-dimensional  point  z  directly  from  the 
joint  MGF  //Z(X)  without  having  to  calculate  or  deal  with  joint 
PDF  pz(z). 
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(27) 


Example  3.  gz(u1, . . . ,uM)  =  U(u1-z1)  •••  U(uM~zM)  , 

where  unit  step  function  U(x)  is  1  for  x  >  0  and  is  0  for  x  <  0. 
Then,  average  (14)  becomes 

00  00 

a  -  I  dUl  •••  J  dUM  PZ(U1 . V 

Z1  ZM 

which  is  the  joint  EDF  of  RV  z  at  argument  z  =  [z^  •••  zMJ' .  The 
gamma  function  (16)  corresponding  to  nonlinearity  (27)  is 


“  e2(  zi '  *  *  *  •  ziq)  •  (  28  ) 


provided  that  the  m-th  Bromwich  contour  Cm  stays  to  the  right  of 
the  origin  for  m=l:M,  and  also  stays  within  the  ROA  of  joint  MGF 
//Z(X).  Relation  (30)  provides  a  method  for  obtaining  the  joint 
EDF  ez(z)  at  arbitrary  M-dimensional  point  z  directly  from  the 
joint  MGF  >t/2(  X)  without  having  to  calculate  or  deal  with  joint 
PDF  p2(z). 
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Example  4.  9Z  (ui  •  •  *  *  'UM)  =  5(u^-z^)  U(Z2~U2)  •••  u^zm-um^  * 
Average  (14)  is  given  by 


_co 


pz(zl'u2 . V  ' 


(32) 


which  is  1  part  PDF  and  M-l  parts  CDF.  M-dimensional  field  point 
z  =  ( z ^  •••  is  arbitrary.  The  gamma  function  (16), 

corresponding  to  nonlinearity  (31),  is 


z2  ZM 

YZ(X)  =  expt-Xj^  z^^ )  J  du2  •••  J  duM  exp(-X2  u2~  •••  -XM  uM)  = 


=  exp(-X'z) 


r  M 

n<-x”) 


if  Re(X  )  <  0  for  m=2:M  . 
m 


(33) 


There  is  no  restriction  on  X1  in  this  particular  gamma  function. 
The  alternative  form  for  average  a  in  equation  (32)  is  given  by 


(  i  2  ii ) 


M 


J 


dX 


A/Z(X) 


exp(-X'z) 


ii'-v 


(34) 


provided  that  the  m-th  Bromwich  contour  Cm  stays  to  the  left  of 

the  origin  for  m=2:M,  and  also  stays  within  the  ROA  of  joint  MGF 

//  (X).  The  particular  contour  C.  for  m=l  can  pass  either  to  the 
z  x 

right  or  left  of  the  origin  in  the  complex  X^  plane  but  must  stay 

within  the  ROA(//  ). 

z 
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(35) 


Example  5.  g  (u, 
2  1 


-  U(um-Z»b 


where  zma  <  zmh  for  m-l:M.  The  difference  of  unit  step  functions 
is  a  unit  pulse  over  interval  (Zma/Zmb).  Average  (14)  is  now 


2lb  zMb 

a  =  J  du^  • • • 

z, 
la 

(36) 

which  is  recognized  as  a  set  of  joint  interval  probabilities. 
Limits  {zmQ}  and  {zmbJ  are  arbitrary.  The  gamma  function  (16) 
corresponding  to  nonlinearity  (35)  is  given  by 


J  duM  Pz(ul- 


'V  - 


Pr ( zla<zl<zlb' 


'WVW 


Ma 


f  z 


M 

rz(X)  =  TT 

m=l 


mb 

J  du  exp(-Xmu) 


,  z 

k  ma 


M  f exp ( -X  z  )  -  exp(-X  z  ,  ) 
T^T  J _ ra  ma _ _  m  mb 7 

=ll  X- 


m 


m 


) 


(37) 


for  all  X.  The  use  of  equation  (37)  in  general  relation  (18) 
yields  the  alternative  average  form 


a 


’  [ dX U{ 


-  ■exP<-x„ii.h)' 


m 


(38) 


where  M-dimensional  contour  C  is  restricted  only  by  being 

required  to  remain  in  the  ROA  of  joint  MGF  //  (X).  These  interval 

z 

probabilities  (38)  can  be  determined  directly  from  the  joint  MGF 
//2(X)  of  RV  z  without  the  need  to  calculate  or  deal  with  the 
joint  PDF  Pz(z)  of  RV  z,  as  would  be  required  by  attempting  to 
use  equation  (36). 
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v.-l  VM-1 

Example  6.  g(u1#...,uM)  =  u1  ’ *  *  UM  for  um  >  0 

and  zero  elsewhere;  parameters  -v  >  0  for  m=l:M.  Let  { zm> , 
m=l:M,  be  positive  RVs  with  joint  PDF  pz  and  joint  MGF  The 

average  (14)  is  given  by 


.  vl-1 
duM  U1 


v  -1 

“  UM  Pz(Ul' *  *  * ,UM) 


(40) 


which  is  recognized  as  the  (vm~l)  th-order  moments  of  RV  z.  The 
gamma  function  (16),  corresponding  to  nonlinearity  (39),  is 


M 

r(X)  -  TT 

m-1 


Ko 


V  -1 

du  u  exp( -A  u) } 

mm  r  mm1 


m  rr(v  )) 

=  rrl-vM 

m=l  L  i  J 


(41) 


m 


if  Re(Xm)  >  0  for  m=l:M.  The  use  of  equation  (41)  in  general 
relation  (18)  yields  the  alternative  form 


a 


T7lr(V) 

m=l 


(42) 


where  denotes  that  the  M  Bromwich  contours  {Cm)  must  all  pass 
to  the  right  of  the  origin  in  their  respective  \m  planes.  There 
is  no  parameter  vector  z  for  this  particular  probabilistic 
problem  ( 40 ) . 
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Example  7.  gfu^/i^)  -  1  for  |u^+u2|  <  1;  zero  otherwise.  (43) 

This  two-dimensional  nonlinearity  is  nonzero  only  within  a  unit 
diamond  D  centered  at  the  origin.  The  average  (14)  is  given  by 

8  =  JJ  dUl  du2  PZ(ul'u2)  =  pr(z1,z2  e  D)  ,  (44) 

D 

which  is  the  probability  of  RV  z  =  [Zl  z2J'  landing  in  the  unit 
diamond  D  located  at  the  origin.  The  gamma  function  (16), 
corresponding  to  nonlinearity  (43),  is 

y(X)  =  JJ  d^  du2  exp(-X1  Uj  -X2  u2)  = 

D 

cosh(X^)  -  cosh(X2) 

=  4  - -  for  all  X1,X2  .  (45) 

A1  X2 

The  use  of  equation  (45)  in  general  relation  (18)  leads  to 
alternative  form 


a 


4 

( i2re) 2 


J  dx  A,Z(X) 

c 


coshfX^)  -  cosh(X2) 


(46) 


where  contour  C  is  restricted  only  to  be  in  the  ROA  of  joint  MGF 
*/z(X). 


Notice  that  nonlinearity  g(u1,u2)  in  equation  (43)  is  not 
separable  in  the  variables  u^,u2;  this  is  in  contrast  to  the 
earlier  six  examples,  where  function  g(ulf...,uM)  was  separable 
in  all  M  variables  {umJ .  Also,  all  the  g-function  examples  here 
have  been  taken  to  be  positive  real  where  they  are  not  zero. 
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SADDLEPOINT  APPROXIMATION 


The  general  M-dimensional  integral  of  interest  is  given  by 
equations  (17)  and  (18)  in  the  form 

a  -  — — =  [  dX  o  (X)  r( X) - - — Si  f  dx  expIA(X)  ]  ,  (47) 

( i2n)M  J  z  (i2a)M  J 

where  the  function  A(X)  in  equation  (47)  is  defined  as 

A(  X)  =  ln[fJz(X)  r(  X)  ]  =  X2(X)  +  In  y(X)  ,  (48) 

and  where  X„(X)  is  the  joint  CGF  (11)  of  RV  z. 
z 


A  real  SP  of  the  total  integrand  in  equation  (47)  is  located 

where  the  M  partial  derivatives  (PDs)  of  exp[A(X)]  are  all  equal 

/\ 

to  zero.  Equivalently,  a  real  SP  X  is  located  where  the  PDs 


9A(  X) 
9X 

m 


0  for  m=l:M  . 


(49) 


Saddlepoint  equation  (49)  constitutes  M  nonlinear  simultaneous 
equations  in  real  variables  {Xm),  m=l:M.  The  location  of  the 

A  A 

real  SP  X  depends  on  both  XZ(X)  and  y(X).  Also,  SP  X  must  lie  in 
the  intersection  ROA(a/„,y)  of  their  two  ROAs.  The  M-dimensional 
contour  C  in  equation  (47)  can  be  taken  to  pass  through  this  real 
SP;  that  is,  each  Bromwich  contour  passes  through  real 

A  A 

component  Xm  of  real  vector  SP  X.  At  this  stage,  there  are  no 
approximations  involved  in  expression  (47);  this  is  an  exact 
result  for  the  average  a  of  interest. 
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DERIVATION  OF  SADDLEPOINT  APPROXIMATION 


Once  the  real  SP  X  in  the  ROA(//z,y)  has  been  located,  expand 
the  A( X)  function  in  equation  (48)  about  that  point  according  to 
M-dimensional  Taylor  series 


M 


A(x>  =  a(  x)  +  y~: 


3A(X) 


__i  3X 
m=l  in 


(Xm  ‘  + 


4  r  A.(m,m)  [\  -  X J  fxm  -  X  J 
2  niT^l  2  "  ^  ra  “  W 


(50) 


where  X^  is  real  and  second-order  coefficients 


,  _v  _  3  A(X) 
A0  (m,m)  =  -r^r  r-T- 


‘2 "  3X  3X 
in  m 


for  m,m=l:M  . 


(51) 


Also,  define  the  MxM  Hessian  matrix  of  A(X)  at  the  SP  X  as 


X2  s  [A2(m,m)]  ,  m,m=l:M  . 


(52) 


Now,  truncate  expansion  (50)  at  second  order,  use  equation 
(49),  and  substitute  the  result  into  equation  (47)  to  get  the 
approximation 


a°  *  Tlfe*  J dX  exp[A,X)  +  * 

At  this  point,  make  the  change  of  variable 

X=  X+  i  for  m=l : M  , 

in  m  m  9 


(54) 


and  define  Mxl  real  vector  s  =  [  s^^  •••  sM]'.  Then,  equation  (53) 
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yields 


a. 


eXptA(M}  1  J  ds  exp[-  |  E  A2(m,m)  s  s  ]  - 

(2n)n  i  L  z  m,m=l  ^  m  -J 


exp[ A( A) ] 


( 2 n ) M/2  [det(A2)] 


*5  ' 


(55) 


where  equation  (52)  has  been  used.  This  result,  a^,  in  equation 
(55)  is  denoted  as  the  zeroth-order  saddlepoint  approximation  to 
exact  result  a.  Observe  that  ag  must  always  be  positive,  since 
MxM  Hessian  matrix  A2  is  always  positive  definite  at  real  SP  A  if 
g(u)  in  equation  (18)  is  real  and  non-negative  for  all  u. 


CORRECTIONS  TO  SADDLEPOINT  APPROXIMATION 

The  first-order  correction  to  ap  is  given  in  reference  1, 
pages  15  -  16,  as 

al  s  aQ  [1  +  cfc]  ,  (56) 

where  c^  is  a  correction  term.  Alternatively,  the  exponential 
version  of  the  saddlepoint  correction  is  (reference  3,  page  180) 

ae  s  ao  exPlct)  »  (57) 

which  has  the  same  first-order  term  as  equation  (56).  However, 
it  has  been  found  numerically  that  approximation  ag  generally 
gives  more  accurate  results  than  either  aQ  or  a^;  accordingly,  ag 
has  been  used  for  the  majority  of  the  numerical  results  to  be 
presented  here. 
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The  correction  term  cfc  uses  the  third-  and  fourth-order 
PDs  of  the  function  A(X),  defined  in  equation  (48);  however, 
these  particular  PDs  only  need  to  be  evaluated  at  the  real 

A 

saddlepoint  X. 

Expressions  were  given  earlier  for  the  gamma  function  y(X) 

corresponding  to  the  joint  CDF,  EDF ,  PDF,  moments,  and  interval 

probabilities.  Also,  the  joint  MGF  //  ( X )  can  be  found  in  closed 

z 

form  for  a  number  of  statistics,  such  as  the  joint  Gaussian, 

quadratic  and  linear  forms  in  dependent  Gaussian  RVs,  and 

linearly  transformed  independent  RVs  with  arbitrary  statistics. 

By  combining  these  results  with  equations  (55)  -  (57),  SPAs  can 

now  be  obtained  for  a  number  of  M-dimensional  probabilistic 

problems  that  were  previously  considered  unsolvable.  These  SPAs 

require  that  the  M-diraensional  field  point,  z,  of  interest  be 

specified  and  that  the  real  SP  X  in  the  ROA(//  ,  y)  be  obtained. 

z 

Strictly  speaking,  closed  forms  for  joint  MGF  fj  (X)  or  joint 

z 

CGF  XZ(X)  are  not  mandatory;  however,  it  must  be  possible  to 

numerically  evaluate  //_ ( X )  or  X„(X)  at  any  M-dimensional  real 

z  z 

point  X  required.  Then,  a  numerical  search  for  the  minimum  of 
total  integrand  /«Z(X)  y(X)  or  XZ(X)  +  In  y(X)  of  equations  (47) 
and  (48),  where  the  search  can  be  confined  to  the  real  {X  }  axes, 

A 

will  locate  the  real  SP  X  in  the  ROA.  The  Hessian  matrix  A2  in 
equation  (52)  will  then  require  second-order  differences  of  A(X) 

A 

be  taken  in  the  neighborhood  of  the  SP  X.  Numerical  evaluation 
of  third-  and  fourth-order  PDs  may  be  questionable,  however. 
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CUMULATIVE  DISTRIBUTION  FUNCTION  OF  THE 
MAXIMUM  OF  M  DEPENDENT  RANDOM  VARIABLES 


RV  z  =  [  z^  •••  Zjj]'  has  statistically  dependent  components 
{zm},  m=l:M.  The  maximum  RV  of  set  {zm}  is  defined  as 

y  =  max(zlf . . . ,zM)  ,  (58) 

which  is  a  random  scalar.  The  first-order  CDF  of  RV  y  is 


cy(y)  =  Pr(y  <  y)  =  Pr(z1  <  y,...,  zM  <  y)  =  cz(y,...,y)  = 


1 

(i2n)M 


J  ax  ,,(X)  exp (-  y  C  \J/TT<-V  > 
r  m=i  /  m=i 

L1 


(59) 


where  M-dimensional  contour  denotes  that  the  Bromwich  contours 
must  all  pass  to  the  left  of  their  respective  origins.  Equation 
(59)  is  a  single  M-dimensional  integral  for  the  CDF  of  maximum  RV 
y  that  allows  for  arbitrary  statistical  dependencies  between  RVs 
{zmJ,  as  reflected  through  the  joint  MGF  /jz(\)  . 


The  A(X)  function  of  equation  (48)  is  given  here  by 


M 


M 


A<X)  =  XZ(X)  -  jC  xm  -  5E  ln<-A_) 


m=l 


m 


m=l 


m 


The  first-order  PDs  of  A(X)  are 


5A(X)  3Xz,X)  1  .  .. 

-sir  ‘  — - y  -  x-  fot  m-l!M 


m 


m 


m 


(60) 


(61) 


The  SP  equations  are  obtained  by  setting  these  M  PDs  to  zero  and 

a 

numerically  solving  for  M-dimensional  SP  X.  However,  the  search 
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for  the  real  SP  solution  of  equation  (61)  must  be  conducted  only 

in  the  left  half  of  each  real  Xm  axis  and  within  the  ROA  of  joint 

MGF  //  (X).  The  search  procedure  must  not  be  allowed  to  wander 

into  any  of  the  right  half  real  Xm  axes;  otherwise,  the  search 

may  drift  off  to  infinity  or  find  a  spurious  real  point  of  the 

analytic  continuation  of  //(X)  outside  the  ROA  (//  ,  y)  that  also 

z  z 

satisfies  the  SP  equations. 


The  second-order  PDs  of  A(X)  follow  from  equation  (61)  as 


32A(X) 

3  X  3  X 
m  m 


9  XZ(X) 
3Xm  3Xm 


_1  , 

.2  mm 
m 


for  m,m=l:M 


(62) 


The  addition  of  the  positive  diagonal  terms  {1/X2}  (at  the  real 
SP)  to  the  Hessian  matrix  of  XZ(X)  serves  only  to  improve  the 
positive  definitenes's  of  the  Hessian  matrix  (62)  for  A(X). 
Equations  (60)  through  (62)  enable  calculation  of  the  zeroth- 
order  SPA  ap  given  by  equation  (55).  The  correction  term  c^.  in 
equations  (56)  and  (57)  also  requires  the  third-  and  fourth-order 
PDs  of  A( X)  and  can  be  derived  from  equation  (62). 


TWO  TEST  EXAMPLES 

Two  test  examples  will  be  used  extensively  in  the  following 
developments.  They  are 

1.  Joint  Gaussian  RVs  {zffi},  m=l:M,  with  arbitrary  Mxl  mean 
vector  and  MxM  covariance  matrix,  and 
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2.  Linearly  transformed  exponential  RVs  {enJ,  n=l:N.  In 

/ 

particular,  RV  z  =  A  e,  where  Nxl  RV  e  is  composed  of 
independent,  identically  distributed  RVs,  each  with  PDF 
exp(-u)  for  u  >  0.  Arbitrary  matrix  A  is  MxN,  thereby  making 
z  an  Mxl  RV  with  dependent  non-Gaussian  components  {zmJ. 

The  following  numerical  examples  all  employ  M  *  4  and  N  -  7. 
Since  there  are  no  analytic  results  for  the  statistics  of  maximum 
variate  y  defined  in  equation  (58),  when  components  {z^}  are 
statistically  dependent,  three  different  methods  are  employed  to 
determine  the  "truth"  for  comparison  with  the  SPAs  that  are 
developed  here;  they  are  Gauss-Hermite  quadrature  for  the  lower 
tail  of  RV  y,  simulation  using  le8  trials  for  the  central  region, 
and  an  asymptotic  form  for  the  upper  tail  of  RV  y. 

To  employ  Gauss-Hermite  quadrature  on  the  exact  M-dimensional 
integral  (47),  the  contour  C  is  first  moved  to  the  real  SP  of 
total  integrand  exp[A(X)].  Then,  an  M-dimensional  linear 
transformation  of  variables  is  utilized  that  makes  the  form  of 
the  new  integrand,  in  the  neighborhood  of  the  SP,  behave  as 

exp[-  |(s \  +  •••  +  4)]  '  (63) 

where  s  is  the  vertical  deviation  from  the  real  SP  in  the 
m 

complex  Xm  plane.  Finally,  the  standard  one-dimensional  Gauss- 
Hermite  relations  are  replicated  in  all  M  dimensions,  using  the 
weighting  factor  (63).  The  number  of  samples  per  dimension  is 
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increased  until  stability  in  the  estimated  integral  is  realized, 
or  until  storage  and  execution  time  problems  become  untenable. 

It  has  been  found  for  the  current  two  test  cases  that  this 
Gauss-Hermite  procedure  is  accurate  enough  only  on  the  lower  tail 
of  maximum  RV  y. 

The  asymptotic  EDF  of  maximum  RV  y  is  given  as  the  sum  of  the 
first-order  EDFs  of  each  of  the  individual  RVs  {zm)  in  equation 
(58).  This  asymptotic  result  becomes  relatively  more  accurate  as 
the  upper  tail  of  RV  y  is  examined.  The  combination  of  the  three 
approaches  above  is  sufficient  to  enable  determination  of  the 
relative  accuracies  of  the  SPAs  over  essentially  the  full  range 
of  arguments  of  the  various  statistics  considered  here. 


GRAPHICAL  RESULTS 

Figure  1  displays  results  for  the  CDF  of  maximum  RV  y  for  the 
Gaussian  test  example  (GTE)  cited  earlier.  The  black  curve  is 
the  simulation  result  (SIM)  using  le8  trials,  while  the  x  points 
are  the  SPAs  obtained  from  exponential  version  ae  in  equation 
(57).  Equations  (59)  through  (62)  were  used  for  these 
computations.  The  red  curve  is  the  asymptotic  result  (ASY)  for 
the  CDF  of  y;  the  asymptotic  result  corroborates  the  simulation 
in  the  central  region  and  significantly  extends  the  CDF 
simulation  on  the  upper  tail.  On  the  lower  tail,  the  asymptotic 
result  is  useless. 
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Figure  2  contains  plots  (for  the  GTE)  of  the  ratio  of  SPA  ag 

to  the  three  attempts  at  "truth",  to  be  called  accuracy  ratios. 

Gauss-Hermi te  quadrature  is  denoted  by  GHQ,  the  simulation  by 

SIM,  and  asymptotic  by  ASY.  The  GHQ  curve  can  be  trusted  for 

threshold  argument  y  less  than  1,  the  SIM  curve  for  y  larger  than 

0,  and  the  ASY  curve  for  y  larger  than  6.  The  end  result  is  that 

the  SPA  a  varies  from  +3%  to  -8%  error  over  the  complete  range 
e 

of  its  argument  values.  The  notation  NG  denotes  regions  where 
the  plotted  results  are  no  good  and  cannot  be  used. 

Figure  3  displays  results  for  the  CDF  of  maximum  RV  y  for  the 
exponential  test  example  (ETE)  cited  earlier.  The  superposed 
points  labeled  with  0  are  the  SPAs  obtained  from  the  zeroth-order 
SPA  ag  in  equation  (55).  Approximation  ag  is  significantly 
poorer  than  exponential  SPA  ag;  in  fact,  ag  becomes  larger  than  1 
for  threshold  argument  y  larger  than  7;  although  SPA  ag  must  be 
positive,  it  need  not  stay  below  1.  By  contrast,  SPA  ag  stays 
below  1  for  all  y  for  these  examples.  An  asymptotic  result  for 
the  CDF  of  RV  y  was  not  computed  for  this  ETE. 

Figure  4  contains  the  accuracy  ratios  for  the  ETE.  The 
zeroth-order  approximation  ag  has  error  variations  from  -70%  to 
+35%,  whereas  ag  only  errs  by  -20%  to  +4%.  This  example 
illustrates  the  need  to  compute  the  correction  term  c^  for  use  in 
equation  (57).  Of  course,  the  additional  effort  required  to 
evaluate  the  third-  and  fourth-order  PDs  needed  for  c^  can  be 
considerable;  see  reference  1  for  an  example. 
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NEED  TO  RESORT  TO  EDF  OF  MAXIMUM  VARIATE  y 


For  scalar  RV  y,  the  EDF  is  simply  ey(y)  =  1  -  cy(y),  from 
which  the  detection  probability  Pd  and  false  alarm  probability  Pf 
of  3  th r e shold— compa r i son  processor  can  be  immediately  calculated 
according  to  Pfl  =  1  -  cy( y | H1 )  and  Pf  =  1  -  cy(y|HQ),  where  H1 
and  H0  are  the  hypotheses  that  a  signal  is  present  and  absent, 
respectively,  if  the  CDF  values  cy(y)  were  exact,  these 
relations  could  be  used  directly.  However,  small  approximation 
errors  in  calculating  the  CDF  values  cy(y)  can  sometimes  result 
in  large  EDF  errors,  depending  on  the  exact  range  under 
consideration. 

For  example,  if  the  exact  Pd  is  0.99,  then  c  (yl^)  =  0.01. 

However,  the  SPA  to  this  CDF  may  be  0.01  (1  +  e),  where  e  is  in 

the  range  +0.1.  This  result  leads  to  an  SPA  for  P.  of  value 

a 

1  -  0.01  (1  +  e)  =  0.99  -  0.01  e  =  (0.989  to  0.991).  This  range 
of  variation  is  probably  quite  acceptable  for  Pd  evaluation. 

On  the  other  hand,  suppose  the  exact  Pf  is  0.001,  meaning 

cy(ylH0)  *  0.999.  But,  for  cy  SPA  0.999  (1  +  e),  the  SPA  to  Pf 

is  1  -  0.999  (1  +  e)  «  0.001  -  0.999  e  =  (-0.099  to  0.101).  This 
range  is  totally  unacceptable.  Thus,  it  is  necessary  to  have  a 

method  for  direct  calculation  of  the  EDF  ey(y)  itself  of  maximum 

RV  y,  especially  for  very  small  EDF  or  Pf  values. 
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EXCEEDANCE  DISTRIBUTION  FUNCTION  OF  THE 
MAXIMUM  OF  M  DEPENDENT  RANDOM  VARIABLES 

The  first-order  EDF  of  maximum  variate  y  in  equation  (58)  is 
not  as  simple  as  that  given  for  the  CDF  in  equation  (59).  The 
added  complexity  for  the  EDF  is  best  illustrated  by  considering 
the  special  case  of  M  =  2;  that  is,  y  =  max(z)  =  max(z1,z2). 

EXCEEDANCE  DISTRIBUTION  FUNCTION  OF  y  FOR  M  -  2 

For  M  =  2,  there  follows  y  =  max(zlfz2).  Consider  a  two- 
dimensional  z^,z2  plane  with  horizontal  and  vertical  lines  drawn 
at  threshold  y  in  both  dimensions.  Label  the  quadrants  centered 
at  the  point  (y,y)  in  the  z^,z2  plane  as  Q^,  Q2,  Q^,  and  Q4,  in 
standard  order.  Then,  CDF  Cy(y)  =  Prfz^  <  y,  z2  <  =  CZ(Y'Y) 

is  the  probability  of  RV  pair  z^,z2  landing  in  quadrant  Q^. 

On  the  other  hand,  EDF  ey(y)  is  the  probability  of  RV  pair 
z^,z2  landing  in  quadrants  Q^,  Q2,  or  Q4 ;  call  these  individual 
quadrant  probabilities  P^  P2,  and  P4,  respectively.  Thus,  there 
follows,  for  the  composite  probability, 

ey(Y>  =  P1  +  p2  +  P4  = 

=  (F1  +  P2)  +  (Pj  +  P4)  -  Px  = 

-  (Px  +  P4)  +  P2  .  (64) 
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The  first  form  in  equation  (64)  adds  up  positive  numbers; 
however,  there  are  three  probability  terms  that  have  to  be 
computed.  The  second  form  in  equation  (64)  (called  inclusion  and 
exclusion)  also  requires  that  three  different  probabilities  be 
calculated;  additionally,  it  requires  the  subtraction  of  a 
positive  quantity.  Finally,  the  third  form  in  equation  (64) 
requires  that  only  two  terms  be  evaluated,  both  of  which  are 
positive.  The  third  form  in  equation  (64)  can  be  interpreted  as 
follows;  >  y)  and  P2  =  Pr^  <  y,  z2  >  y).  The 

P2  term  is  a  mixed  probability,  namely,  1  part  CDF  and  1  part 
EOF. 

If  the  quadrant  probabilities  { Pfc } ,  k=l:4,  are  not  calculated 
exactly,  but  perhaps  obtained  by  SPAs,  the  second  form  in 
equation  (64)  will  be  subject  to  possible  loss  of  significance 
due  to  the  negative  term.  Therefore,  this  second  form  is  not 
recommended  for  use  when  only  approximations  to  the  individual 
probabilities  are  available.  Also,  forms  1  and  2  in  equation 
(64)  require  three  probability  evaluations,  whereas  form  3 
requires  only  two  probability  evaluations;  this  will  become 
significant  in  M  dimensions  for  M  >  2. 
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EXCEEDANCE  DISTRIBUTION  FUNCTION  OF  y  FOR  M  >  2 

For  general  M,  the  maximum  variate  is  given  by  equation  (58) 
as  y  =  max(Zj,...,2M).  The  EDF  of  RV  y  can  be  written  as  an 
obvious  generalization  of  the  third  form  of  equation  (64): 

ey(y)  =  Pr(y  >  y)  =  Pr^  >  y)  +  Pr(z1  <  y,  z2  >  y)  + 

+  Pr(z^  <  y,  z2  <  y,  z^  >  y)  +  ••• 

+  Pr(zx  <  y,,..,  zM_1  <  y,  zM  >  y)  .  (65) 

Equation  (65)  contains  M  terms,  all  of  which  are  positive;  there 
are  no  cancellations  involved  in  this  form.  Whereas  the  number 
of  terms  in  form  (65)  increases  only  linearly  with  dimensionality 
M,  by  contrast  the  corresponding  generalizations  of  forms  1  and  2 
in  equation  (64)  would  involve  2M  -  1  terms.  Thus,  the  effort 
here  will  employ  form  (65)  for  all  future  calculations. 


By  reference  to  equations  (5),  (8),  and  (30),  the  EDF  in 
equation  (65)  can  be  written  as  the  sum  of  a  number  of  contour 
integrals: 


Vy)  =55fakl 

cr 
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)  exp(-y\1) 
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where  auxiliary  MGFs  are  defined  at  stage  n  as 


(  X^ '  *  *  *  '  )  —  ^  z  ^  ^1 '  *  *  *  '  ^n '  ^  ^ '  ’  *  *  '  ^  ^  for  n=l;M  .  (67) 

Stage  variable  n  varies  from  1  to  M.  There  immediately  follows, 
from  equation  (67),  the  analogous  auxiliary  CGFs: 

Xjj  (  X^  ,  •  •  <  ,  (  \  i  •  •  •  /  X^  ,0,0, ...,0)  for  n=l;M  •  (68) 

The  result  in  equation  (66)  is  exact.  In  order  to  calculate 
the  EDF  ey(y)  via  this  equation,  it  is  necessary  to  conduct  one 
one-dimensional  integral  plus  one  two-dimensional  integral  plus 
...  plus  one  M-diraensional  integral.  To  achieve  the  SPAs,  this 
will  entail  determination  of  a  one-dimensional  real  SP  plus  a 
two-dimensional  real  SP  plus  ...  plus  an  M-dimensional  real  SP. 
Alternatively,  the  first  few  integrals  in  equation  (66)  could  be 
conducted  rather  accurately  via  FFTs ,  if  desired. 

Since  all  the  probabilities  in  equation  (66)  are  positive, 
the  corresponding  SPAs  are  also  positive.  When  a  number  of 
positive  approximations  are  added,  the  general  effect  is  to 
average  the  errors  of  each,  and  to  realize  some  improvement  of 
the  accuracy.  An  exception  occurs  when  one  of  the  approximations 
has  a  very  large  error,  in  which  case  it  may  dominate  the  overall 
error  of  the  sum.  The  least  accurate  SPAs  to  the  individual 
terms  in  equation  (66)  are  not  known.  The  approach  utilized  here 
was  to  evaluate  the  SPAs  for  all  the  terms  in  equation  (66)  and 
add  them;  no  FFTs  were  employed. 
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At  the  n-th  stage  of  equation  (66),  the  A(X)  function  of 


equation  (48)  takes  the  form 


An(Xl' ‘ • *'Xn)  "  Xn(Xl' 


*  * '  Xn  > 


n  n-1 

y  C  xm  -  ZZ  in<-V 


m=l 


m=l 


ln(Xn)  . 


(69) 


The  PDs  are  given  by 
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for  m=l:n 


(70) 
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Therefore,  the  nxl  SP  vector  X^n^  =  [xjn^***  X^n^J'  at  the  n-th 
stage  is  the  solution  of  the  n  simultaneous  nonlinear  equations 


a  Xjj  (  r ...  r  xn ) 


ax_ 


m 


(n) 


0  for  m=l:n  . 


(71) 


Since  the  last  contour  in  each  integral  in  equation  (66)  must 
pass  to  the  right  of  the  origin  in  the  complex  Xn  plane,  the  n-th 
solution  component,  X^  '  in  equation  (71),  must  be  positive  real. 
Conversely,  all  the  other  n-1  components  of  the  SP  vector  X^n^ 
must  be  negative  real. 


In  order  to  effect  the  calculation  of  equation  (66)  for  the 
total  EDF  ey(y),  it  is  necessary  to  solve  equation  (71) 
repeatedly,  as  stage  number  n  varies  from  1  to  M.  Thus,  the 
dimensionality  of  the  SP  search  varies  from  n  «  1  dimension  to 
n  «  M  dimensions.  Then,  all  M  SPAs  to  all  the  terms  in  equation 
(66)  are  added  to  yield  the  final  approximation  to  the  EDF  ey(y). 
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GRAPHICAL  RESULTS 


Figure  5  displays  the  EDF  ey(y)  of  maximum  RV  y  for  the  GTE, 
as  determined  from  equation  (66).  The  curves  have  the  same 
identifications  as  the  CDF  curves  did  earlier.  The  results  for 
the  asymptotic  behavior  and  the  SPA  ag  overlie  each  other  well 
into  the  upper  tail,  while  the  simulation  estimates  eventually 
become  unstable  due  to  an  insufficient  number  of  trials  (le8). 

Figure  6  gives  the  corresponding  accuracy  ratios  for  this 
GTE.  The  SPA  is  off  by  -3%  on  the  lower  tail;  this  error 
increases  to  +1%  for  threshold  y  near  1;  then,  the  error  drifts 
down  to  almost  -3%  before  turning  up  for  large  y  arguments. 

There  is  an  imprecise  transition  region  for  y  in  the  range  (8,13) 
where  the  simulation  results  are  gradually  replaced  by  the 
asymptotic  results.  The  NG  regions  indicate  where  the  curves 
should  definitely  be  ignored.  Again,  Gauss-Hermite  quadrature, 
simulation,  and  asymptotic  behavior  were  employed  to  cover  the 
entire  range  of  argument  values.  The  small  errors  for  this  GTE, 
namely  +3%,  are  very  encouraging;  this  may  be  partially  due  to 
the  averaging  effects  of  the  individual  SPAs  mentioned  above. 

Figure  7  displays  the  EDF  results  for  the  ETE.  The  SPA  a 

e 

and  the  asymptotic  results  overlie  each  other  well  into  the  upper 
tail.  The  accuracy  ratios  in  figure  8  vary  from  -8%  near  y  =  0, 
to  +6%  for  y  near  3.  The  simulation  results  in  figures  7  and  8 
for  larger  y,  that  is,  the  upper  tail,  cannot  be  trusted,  due  to 
instability.  No  asymptotic  results  were  obtained  for  this  ETE. 
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PROBABILITY  DENSITY  FUNCTION  OF  THE 
MAXIMUM  OF  M  DEPENDENT  RANDOM  VARIABLES 

The  interest  in  this  section  is  on  the  first-order  PDF  of  the 
maximum  variate  y  =  max( z^ , . . . , zM) .  One  obvious  possibility  is 
to  take  the  derivative  of  the  first-order  CDF  Cy(y)  in  equation 
(59)  with  respect  to  y.  The  result  is 


py(y) 


( i2n) 
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dX 


"z(X) 


.  M  .  M  /  M 

exp  -yC  \J  C  (-Xm)/rT(-Xm 

m=  1  m=  l  /  m=  I 


),  (72) 


where  the  contours  C^  must  pass  to  the  left  of  the  origins  in 
each  complex  Xm  plane. 


Direct  use  of  equation  (72)  yielded  a  rather  poor  SPA  to  the 
PDF  of  y.  This  is  believed  to  be  due  to  the  factor  involving  the 
sum  of  (-Xm)  terms  in  the  numerator  of  the  integrand  of  equation 
(72).  The  joint  MGF  fj (X),  defined  in  equation  (9),  can  be  seen 
to  decay  in  magnitude  as  X  deviates  from  the  M-dimensional  real 
SP  X  on  a  Bromwich  contour;  the  SP  is  located  on  the  real  axes  of 
the  complex  planes  {X  i,  essentially  because  PDF  p,  is  a  positive 
real  function  where  it  is  nonzero.  Likewise,  the  product  in  the 
denominator  of  equation  (72)  increases  in  magnitude  as  X  deviates 
from  the  real  SP.  Meanwhile,  the  exp  factor  in  equation  (72) 
maintains  constant  magnitude  as  X  deviates  from  the  real  SP  along 
a  Bromwich  contour.  The  effect  of  the  combination  of  these  three 
factors  alone  would  be  to  cause  the  magnitude  of  the  integrand  of 
equation  (72)  to  decrease  as  X  deviates  from  the  real  SP;  this  is 
a  favorable  situation  in  which  to  develop  an  SPA. 
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However,  the  sum  of  the  (-Xm)  terms  in  the  numerator  of 
equation  (72)  increases  as  X  deviates  from  the  real  SP  along  a 
Bromwich  contour.  This  increase  in  magnitude  counters  the 
desired  decay  of  the  total  integrand  in  equation  (72),  which  is 
required  in  order  to  achieve  a  reasonably  accurate  SPA.  A 
numerical  example  employing  equation  (72)  will  be  displayed 
shortly. 


An  approach  that  eliminates  the  troublesome  summation  factor 
from  the  numerator  of  equation  (72)  is  to  cancel  each  of  the  M 
terms,  -Xm,  with  the  corresponding  term  in  the  denominator 
product.  The  net  result  can  then  be  written  as 


M 

pv(y)  =  YZ 

y  n=l 


where  all  the  integrals  at  the  n-th  stage  are  M-dimensional . 

Also,  the  M-dimensional  contour  changes  with  stage  number  n. 

Component  contours  C £n)  pass  to  the  left  of  the  origin  for  m=l:M, 

m-p^n,  while  contour  is  arbitrary.  of  course,  must  also 

stay  within  the  ROA(*/  )  for  each  n=l:M. 

z 


The  major  difficulty  with  equation  (73)  is  that  it  requires  M 
M-dimensional  integrals,  whereas  equation  (72)  only  required  one 
M-dimensional  integral.  Thus,  it  will  be  necessary  to  locate  M 
M-dimensional  real  SPs  for  equation  (73).  Also,  the  question 
arises  as  to  whether  all  the  terms  in  equation  (73)  are  positive. 
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DIRECT  PROBABILITY  DENSITY  FUNCTION  DETERMINATION 


A  direct  determination  of  the  PDF  of  maximum  RV  y,  which  is 
guaranteed  to  have  all  positive  terms,  will  now  be  derived. 
Consider  threshold  value  y  and  an  infinitesimal  increment  dy 
located  at  y.  Then,  if  py(y)  is  the  PDF  of  maximum  RV  y,  the 
quantity  py(y)  dy  can  be  interpreted  as  the  probability  that  the 
largest  RV  y  lies  in  the  interval  (y,y+dy).  But,  there  are  M 
ways  that  this  event  can  occur:  RV  zm  could  lie  in  this  interval 
while  all  M-l  of  the  other  RVs  {zm}  lie  below  threshold  y.  Since 
m  can  vary  from  1  to  M,  it  is  necessary  to  sum  over  the 
probabilities  of  each  event,  getting  total 

Py(y)  dy  =  Pr ( y  <  y  <  y+dy)  = 

=  Pr ( y  <  z1  <  y+dy,  *2  <  y, . . . ,  *M  <  y)  +  • • • 

+  Pr^  <  y,...,  z^  <  y,  y  <  zM  <  y+dy)  .  (74) 

All  M  terms  in  equation  (74)  are  certainly  positive,  being 

probabilities.  The  probabilities  of  events  involving  more  than 

2 

one  hit  in  the  interval  of  width  dy  are  of  order  dy  or  higher, 
and  are  therefore  negligible. 

To  illustrate  the  transformation  of  equation  (74)  into  the  X 
domain,  consider  the  first  term  and  define  nonlinear  function 

g1(u1, .  . .  ,uM)  =  [Ufu^y)  -  U(u1-y-dy)]  U(y-u2)  •••  U(y-uM)  .  (75) 
Then,  the  average 


39 


y+dy  y 

al  -  J  du  Pz(u)  g1(u)  -  J  dux  J  du2 

y  —oo 


y 

J  dUM  Pz^Ul'***UM^ 


is  immediately  recognized  as  the  first  term  on  the  right-hand 
side  of  equation  (74).  The  gamma  function  corresponding  to 
nonlinear  function  g1  in  equation  (75)  is 


y+dy  y 

T1(X)  *  J  dul  J  du2 

y  —co 


y 

•  • •  J  duM  exp(-X'u)  = 


(77) 


if  Re(Xro)  <  0  for  m=2:M.  Substitution  of  equation  (77)  into 
equation  (18)  yields  the  alternative 


_d£. 


(i2u) 


M 


dX  yuz(X)  exp  (y  } 

JL, 

1) 

i=l  / 

TT(-V 

m=2  m 


(78) 


provided  that  contour  passes  to  the  left  for  m=2:M,  while 

contour  cj1*  is  arbitrary,  but  within  the  ROA(//  ). 

z 


Continuing  this  procedure  for  the  remaining  terms  in  equation 
(74)  leads  to  the  result  Py(y)  dy  **  a^  +  •  •  •  +  a^,  at  which  point 
the  cancellation  of  common  factor  dy  leads  precisely  to  equation 
(73).  Thus,  it  can  be  concluded  that  all  M  integrals  in  equation 
(73)  will  yield  positive  values  and,  therefore,  positive  SPAs. 
Equation  (73)  will  require  a  solution  for  M  M-dimensional  real 
SPs  and  an  evaluation  of  M  M-dimensional  integrals  for  the  SPAs. 
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SADDLEPOINT  EQUATIONS 

At  the  n-th  stage,  n=l:M,  the  A(X)  function  of  equation  (48) 
takes  the  form,  from  equation  (73), 


An(X)  =  XZ(X) 


M 

y  IZI  \ 

m=l 


m 


S  in,'x») 

m=i 
m^n 


(79) 


The  PDs  are  then 


3 A  (  X)  3X_(  X) 

n  z 


ax 


m 


ax 


m 


-  y  -  t-  +  t- 5 


m 


n 


mn 


for  m=l:M  , 


(80) 


leading  to  the  SP  equations  at  the  n-th  stage. 


axz(X) 


ax_ 


m 


(n) 


-  y  -  - 


(n) 


m 


+  ■s-r- — r  S  =  0 
•>  ( n )  mn 

n 


for  m=l:M 


(81) 


The  M-dimensional  real  SP  vector  X^n^  must  have  X^n^  <  0  for 

IU 

^  /  n  \ 

m*l:M,  m^n;  however,  component  X^  ’  is  arbitrary,  but  within  the 
ROA(//z).  Also,  the  solution  of  equation  (81)  must  be  conducted 
repeatedly  for  stage  number  n=l:M.  Thus,  M  M-dimensional  SPs 
must  be  computed  and  M  M-.dimensional  SPAs  obtained;  when  added 
together,  they  will  constitute  the  total  SPA  to  the  PDF  of 
maximum  RV  y. 
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GRAPHICAL  RESULTS 


Figure  9  presents  the  simulation  and  SPA  results  for  the  PDF 
of  maximum  RV  y  for  the  GTE,  in  addition  to  the  asymptotic 
behavior.  Whereas  the  simulation  results  become  unstable  for 
threshold  y  greater  than  10,  the  SPA  ag  and  asymptotic  curves 
continue  to  track  each  other  well  down  on  the  upper  tail.  On  the 
lower  tail,  the  SPA  tracks  the  rapid  drop-off  of  the  PDF  while 
the  simulation  becomes  untrustworthy. 

Figure  10  is  a  repeat  of  figure  9  except  that  the  scales  have 
been  changed  and  the  SPA  of  the  single  integral  approach  in 
equation  (72)  has  been  added,  as  labeled  with  symbol  O.  Although 
this  latter  SPA  is  excellent  on  the  lower  tail,  it  deviates 
significantly  on  the  upper  tail,  yielding  useless  results  in  that 

region.  A  possible  reason  for  this  behavior  was  discussed  in  the 
sequel  to  equation  (72). 

Figure  11  gives  the  accuracy  ratios  for  this  GTE.  By  piecing 
together  the  GHQ,  SIM,  and  ASY  results,  it  is  estimated  that  the 
SPA  for  the  PDF  of  RV  y  is  about  +5%  in  error  over  the  complete 
range  of  threshold  y.  There  is  an  additional  curve  labeled 
ae/dpDF,  which  is  an  attempt  at  approximating  the  PDF  by  using 
local  differences  of  the  corresponding  CDF.  Although  adequate  on 
the  lower  tail,  the  dPDF  approximation  quickly  becomes  unusable 
for  moderate  y  values;  the  explanation  for  this  degradation  will 
be  presented  shortly. 
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Figure  12  displays  the  PDF  results  for  the  ETE.  Again,  the 

SPA  ag  has  no  trouble  tracking  the  PDF  on  both  tails  as  well  as 

the  central  region.  The  simulation  results  utilized  a  bin  width 

of  0.1,  thereby  realizing  relatively  few  hits  per  bin  on  the 

upper  tail,  at  least  for  the  le8  trials  used.  The  SPA  a  can  be 

e 

seen  tracking  right  through  the  middle  of  a  smoothed  version  of 
the  simulation  PDF  estimate. 

Figure  13  contains  the  accuracy  ratios  for  this  ETE.  On  the 
lower  tail,  the  SPA  is  in  error  by  -30%,  while  in  the  central 
region,  the  relative  error  is  about  +5%.  The  simulation 
estimates  on  the  upper  tail  cannot  be  trusted;  no  asymptotic 
results  were  derived  for  this  example. 


DANGER  OF  USING  DIFFERENCES  OF  CDF  SPA  VALUES 


Due  to  the  large  amount  of  effort  required  to  evaluate  the 
PDF  Py(Y)  RV  Y  by  means  of  equation  (73),  an  attempt  was  made 
to  estimate  the  PDF  by  taking  local  differences  of  the  CDF  Cy(y) 
in  equation  (59);  the  latter  formula  only  requires  one 
M-dimensional  integral  instead  of  M  M-dimensional  integrals.  To 
ascertain  the  limitations  of  this  approach,  consider  the  PDF 
approximation  at  argument  x  obtained  by  the  ratio 


c(x+A)-c(x-A)  _  ,  . 

2A -  =  Pa(x) 


(82) 
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If  CDF  c  is  exact,  the  error  in  approximation  (82)  is  of  the 

2  9 

order  of  A  ;  that  is,  Pa(x)  =  p(x)  +  0(A),  where  p(x)  is  the 

exact  value  of  the  PDF  at  x. 


However,  if  the  CDF  values  are  themselves  approximations,  the 
estimated  PDF  takes  the  form 

(1  +  e.)  c ( x  +  A)  -  (1  +  e_)  c(x  -  A) 

- 2A -  s  Pb(x)  •  (83) 

This  can  be  expanded  as 


€2  C  4  4*  C  a 

PbU)  =  pa(x)  +  - 2A —  c{x)  +  - 2 — ~  p(x)  +  0(eA)  *  (84) 


Several  possibilities  exist.  If  si  =  e2  *  0,  then 
Pb(x)  =  pg(x)  given  in  equation  (82).  On  the  other  hand,  if 
el  =  e2  =  c  **  0,  then 

Pb(x)  *  Pa<x>  +  c  P(x>  =  P(x)  (1  +  e)  ,  (85) 

which  is  acceptable.  However,  if  jt  0,  and  CDF  value  c(x) 

is  not  small,  then  the  term  ( -  €2)  c(x)/(2A)  in  equation  (84) 
can  be  large  relative  to  p(x);  observe  the  division  by  the  small 
quantity  A.  in  fact,  on  the  upper  tail  of  an  RV,  the  CDF  c(x)  is 
approaching  1  while  the  PDF  p(x)  is  approaching  zero,  making  the 
situation  progressively  worse.  Thus,  the  use  of  equation  (83) 
for  approximating  the  PDF  will  certainly  develop  problems  on  the 
upper  tail,  and  may  develop  problems  earlier.  Figure  11  is  an 
illustration  of  this  effect. 
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loglO(PDF)  loglO(PDF) 


Figure  9.  PDF  of  Maximum  RV  for  GTE 


Figure  10.  Shortcut  PDF  of  Maximum  RV  for  GTE 


Figure  13.  Ac 


STATISTICS  OF  THE  RANGE  VARIATE 


Suppose  RV  x  =  [xj^  •••  xN]'  is  a  set  of  N  dependent  RVs  with 
joint  MGF  ^x( fOjj) .  The  range  variate  for  this  set  of  N  RVs 
is  defined  as 

y  =  max(x^ , . . . f xN)  -  min(x^, . . . ,3^)  .  (86) 

The  first-order  CDF,  EDF ,  and  PDF  of  random  scalar  y  are  of 
interest. 


MATHEMATICAL  MANIPULATIONS 

In  order  to  illustrate  the  procedure  to  be  employed,  a 
numerical  example  is  used;  namely,  N  is  taken  as  4,  in  which 
case  equation  (86)  becomes 

y  =  max(x^ , . . . ,x4 )  -  min(x^, . . . ,x4)  .  (87) 

Equation  (87)  can  be  manipulated  into  a  familiar  form  by  defining 
the  M  =  12  (=  N(N-l))  difference  variables 

Z1  =  xi  -  *2  '  Z4  =  X2  "  X1  '  Z7  x3  “  X1  '  Z10=  X4  “  X1  ' 

z2  =  xi  -  x3  ,  *5  =  x2  -  x3  ,  z8  =  x3  -  x2  ,  z11=  x4  -  x2  , 

z3  =  X1  "  x4  '  Z6  =  X2  “  X4  '  z9  =  X3  “  x4  '  Z12=  x4  _  x3  ' 

(88) 

There  follows,  from  equations  (87)  and  (88), 

y  =  max( z^ , . . . , zM)  ,  M  =  12  ,  (89) 


49 


while  linear  transformation  (88)  can  be  expressed  as 


z  =  A  x  . 


(90) 


Matrix  A  is  12x4  (MxN)  and  has  rank(A)  =  3  (N-l).  This  can 
be  deduced  from  equation  (88)  by  noting  that  only  three  of  the  12 
RVs  *zm*  can  be  indePendently  specified;  the  remainder  follow  as 
simple  differences.  (The  remaining  unlisted  four  differences 
that  are  possible  in  equation  (87)  are  all  zero  and  can  never  be 
the  range  variate.)  Therefore,  the  joint  PDF  of  RV  z  would 
involve  a  conditional  component  with  9  ((N-l)2)  delta  functions. 
Nevertheless,  SPAs  can  still  be  developed  for  the  statistics  of 
the  range  variate  y,  as  given  by  equations  (89)  and  (90). 


Form  (89)  was  encountered  earlier  in  equation  (58).  The  CDF 
of  this  form  of  RV  y  was  given  in  equation  (59),  along  with  the 
A(X)  function  and  its  low-order  PDs  in  equations  (60)  through 
(62).  For  easy  reference,  this  last  equation  is  repeated  here: 


32A( A) 

3  A  3  A 
m  m 


32Xz(X)  1 

VT  3X  +  "2  Snun  £or  “'S'1'" 

m  m  A  — 


(91) 


Although  the  MxM  matrix  of  second-order  PDs  of  x  (A)  has  rank 

z 

3  (N-l),  with  three  positive  eigenvalues,  the  matrix  of  second- 
order  PDs  of  A( A)  has  the  full  rank  12  (M)  due  to  the  {1/A2} 
diagonal  terms,  and  the  Hessian  matrix  A2  is  positive  definite. 
Thus,  the  determinant  encountered  in  SPA  aQ  in  equation  (55)  is 
positive,  and  the  SPA  is  well  behaved. 
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The  joint  MGF  //^ ( X)  required  for  equation  (91)  can  be 
obtained  upon  use  of  equation  (90),  according  to  expectation 

/t/z(X)  =  E  exp(z'X)  =  E  exp(x'A'X)  =  */  ( A'X)  .  (92) 

The  corresponding  joint  CGF  is  XjX)  =  X„(A'X).  The  vector  X  is 
Mxl  while  vector  A'X  is  Nxl,  which  is  the  dimensionality  of  input 
RV  x  in  equation  (86). 

A  major  drawback  with  form  (89)  for  the  range  variate  y  is 
that  the  original  set  of  N  RVs  in  equation  (86)  has  blossomed  to 
the  larger  number  of  RVs,  M  =  N(N-l).  Therefore,  the  SP 
equations  (61)  now  number  M  instead  of  N,  making  this  a  very 
difficult  numerical  problem  for  large  M.  For  example,  if  N  =  10, 
then  M  =  90,  which  means  that  a  search  in  90-dimensional  X-space 
is  required.  An  attempt  at  reducing  the  number  of  RVs  below  the 
number  M  =  N(N-l)  required  by  formulation  (88)  was  attempted  by 
taking  absolute  values  of  the  differences.  Although  this  reduced 
the  number  of  RVs  {zm}  to  N(N-l)/2,  the  resulting  SPA  to  the  CDF 
of  RV  y  was  poorer  than  the  approach  using  equation  (88);  the 
reason  for  the  degradation  was  not  explored. 

The  form  (89)  for  range  variate  y  was  also  encountered  in  the 
considerations  for  the  first-order  EDF  and  PDF  of  the  maximum 
variate;  therefore,  all  the  earlier  results  for  SP  equations  and 
PDs  of  the  EDF  and  PDF  can  be  brought  to  bear  directly  on  the 
range  variate,  with  due  consideration  again  given  to  the 
limitations  imposed  by  a  large  value  of  M  =  N(N-l). 
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GRAPHICAL  RESULTS 


Figure  14  displays  results  for  the  CDF  of  range  variate  y  in 
equation  (89)  for  the  GTE.  The  SPA  aQ  tracks  the  lower  tail  well 
below  the  capability  of  the  simulation  (le8  trials),  but 
underestimates  the  true  CDF  values  in  the  central  region  and 
upper  tail.  The  EDF  SPA  ag  in  figure  15  is  more  accurate  over 
its  entire  range,  except  at  the  lower  tail.  Finally,  the  PDF  SPA 
ae  in  figure  16  gives  an  accurate  representation  of  both  tails, 
but  is  a  slight  underestimate  in  the  central  region. 

The  corresponding  results  for  the  ETE  are  given  in  figures  17 
through  19.  The  accuracies  of  the  CDF,  EDF,  and  PDF  appear  to  be 
slightly  better  than  those  for  the  GTE  above.  Despite  the  sharp 
drop-off  of  the  PDF  in  figure  19  at  the  lower  limit,  the  SPA 
continues  to  track  the  decay;  the  near-origin  threshold  value  is 
y  =  0.001,  for  which  loglO(PDF)  =  -7.5. 
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Threshold  y 
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Figure  17.  CDF  of  Range  Variate  for  ETE 


loglO(POF)  loglO(EDF) 


Figure  18.  EDF  of  Range  Variate  for  ETE 


Figure  19.  PDF  of  Range  Variate  for  ETE 
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THREE  ADDITIONAL  PROBABILISTIC  PROBLEMS 


The  following  three  problems  illustrate  some  additional 
capabilities  of  the  approach  using  the  joint  MGF  in  the  transform 
domain.  Outlines  of  the  solutions  are  given,  but  no  numerical 
results  have  been  evaluated. 


JOINT  CDF  OF  THE  TWO  LARGEST  RANDOM  VARIABLES 

Random  vector  z  =  [z^  •••  zMlf  is  observed;  these  RVs  are 

statistically  dependent  with  joint  MGF  //  (X).  This  data  vector  z 

z 

is  ordered  into  set  y  =  [y^  •••  yMl'  where 

*M  *  *(1-1  s  •••  *  y2  *  yi  •  <93> 

The  joint  CDF  of  the  two  largest  RVs,  y^  and  y2,  is  of  interest; 
that  is,  probability 

<  y±'  Y2  <  with  thresholds  y2  <  y^^  (94) 

is  desired  for  arbitrary  joint  MGF  jj  (X). 

z 

There  are  two  ways  event  (94)  can  occur.  In  the  first,  all 
the  RVs  { zm }  can  be  less  than  threshold  y2.  This  is  given  by 

Pr(yx  <  y2,  y2  <  y2)  =  cz(y2,y2, . . . ,y2)  ,  (95) 

which  is  simply  the  joint  CDF  of  RV  z  at  common  argument  y2. 

The  second  way  event  (94)  occurs  is  when  the  m-th  RV  z  lies  in 
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the  interval  (y2»Y1)  while  all  the  other  RVs  {zmJ  lie  below 
threshold  y2 .  Since  RV  number  m  can  vary  from  1  to  M,  the  total 
probability  of  this  event  is  given  by  sum 


M 

Pr(y2  <  y1  <  ylf  y2  <  y2)  -  5ZU  Pr(y2  <  zm  K  Yl}  zn  K  y2'  n  ^  m>- 


Reference  to  equation  (18)  reveals  that  nonlinear  function  g(u) 
should  be  selected  as  the  form 


g(u 


1' 


'V 


S  Hum  *  y2)  *  °K  -  yl)j  J1  o(y2  -  uj.  (97) 
m=i  n=i 


n^m 


The  corresponding  gamma  function  is  then 


r(  X- 


'V 


eXp(~Xmy2)  ~  exP(-Vl) 
i=l  Xm 


M  ( 

TT 

n=l  l 
n^m 


exp(-Xny2) 


-X 


(98) 


n 


For  a  given  value  of  m  in  the  outer  sum,  Xm  is  unrestricted  but 
the  parameters  { Xn >  in  the  inner  product  must  satisfy  Re(Xn)  <  0 
for  all  n  ?  m.  Of  course,  all  { Xm)  must  always  remain  within  the 
ROA  (aO  of  the  joint  MGF  /u_(X).  Expression  ( 98 )  can  now  be  used 
in  equation  (18)  to  evaluate  probability  (96). 


To  get  the  CDF  for  the  special  case  of  the  second-largest  RV 

y2,  simply  let  threshold  y1  =  +®  above.  Then,  c  (y2)  is  the 

2 

probability  Pr(y2  <  y 2)*  It  can  be  appreciated  that  the 
probability  of  obtaining  the  second-order  CDF  of  the  third- 
largest  and  seventh-largest  RVs  of  set  z  would  be  extremely 
formidable  due  to  the  excessive  number  of  cases  to  evaluate. 
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OR-ING  AND  POST-AVERAGING 


The  or-ing  operation  consists  of  taking  the  maximum  of  a  set 
of  RVs  and  comparing  with  a  threshold  for  purposes  of  deciding  on 
signal  presence  or  absence.  However,  for  low  signal-to-noise 
ratios,  it  is  frequently  desirable  to  perform  further  averaging 
prior  to  the  threshold  comparison.  This  leads  to  consideration 
of  the  random  quantity 

y  =  a  max ( , . . . , a^)  +  0  max(b^, . . . ,bg)  +  max(c^, . . . ,c^)  + 

(99) 

where  the  RVs  a,  b,  c,...  could  be  statistically  dependent  on 
each  other.  Expression  (99)  can  be  modified  into  the  form 

y  =  max(zlf . . . ,zM)  ,  (100) 

where  M  ■  A  B  C  ...  The  reason  for  this  product  is  because  any 
element  of  the  RV  a  must  be  allowed  to  interact  with  any  element 
of  the  RV  b,  etc.,  in  equation  (99).  The  RV  z  in  equation  (100) 
is  simply  a  linear  transformation  of  the  RVs  a,  b,  c, . . .  and  its 
joint  MGF  can  be  determined  from  the  joint  MGF  of  RVs  a,b,c,...; 
however,  a  note  of  caution  is  in  order  here:  dimension  size 
M  =  A  B  C  ...  can  get  very  large  very  quickly,  being  the  product 
of  the  sizes  of  the  RVs  in  equation  (99). 

The  minimum  of  a  set  of  RVs  can  be  easily  translated  to  more 
familiar  terms  by  using  the  relation 

min( , . . . ,vN)  =  -  max( -v^ , . . . , -vN)  .  (101) 
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RATIO  OF  DEPENDENT  RANDOM  VARIABLES 


Let  RV  z  =  [ z ^  z2l'  have  joint  MGF  pz(  a-^c^)  with  ROA(//z). 

The  statistics  of  the  ratio  of  RVs,  namely, 

w  "  zl/z2  '  (102) 

are  of  interest.  In  this  subsection,  RV  z2  is  presumed  to  be 
positive;  the  general  case  of  arbitrary  z2  values  is  treated  in 
the  appendix. 

The  first-order  CDF  of  RV  w  is  given  by 
cw(w)  =  Vx(z1/z2  <  w)  »  Pr(z1  <  w  z2)  =  Pr^  -  w  z2  <  0)  .  (103) 

Define  auxiliary  RV  v  =  z1  -  w  z2 .  Then, 

cw(w)  =  Pr ( v  <  0)  =  cy(0)  =  jjz  J  dX  //v(X)/(-X)  ,  (104) 

C1 

which  is  a  one-dimensional  contour  integral  passing  to  the  left 
of  the  origin  in  the  complex  X-plane.  The  required  first-order 
MGF  of  RV  v  is 

^V(X)  =  E  exp(  X  v)  «  E  exp(  X  z^  -  X  w  z2 )  =  /«2(X,  -wX)  ,  (105) 

where  (X,  — wX)  e  ROA (//z).  Substitution  in  equation  (104)  yields 
the  first-order  CDF  of  ratio  w  in  the  form 

Cw<w)  -  IK  J  dX  "z(X'  -wX)/(-X)  .  (106) 
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The  corresponding  first-order  EDF  of  ratio  w  can  be  found  in 
a  similar  manner: 

ew(w)  "  Til  J  dX  //z(X'  "wX>/x  '  (107) 

Cr 

from  which  there  follows  the  first-order  PDF 

pw(w)  =  na  J  dX  *42)<x'  "wX)  *  aos) 

c 

The  superscript  (2)  denotes  a  partial  derivative  of  joint  MGF 
/t/z^al,<x2^  w*th  respect  to  the  second  argument  a2.  The  contour  C 
in  integral  (108)  need  satisfy  only  the  constraint  (X,  -wX)  e 
ROA(//z)  whereas  contour  integrals  (106)  and  (107)  have  additional 
half-plane  restrictions.  Additional  details  about  the  A(X) 
function  for  equation  (108)  are  found  in  the  appendix. 

If  RVs  z^  and  Z2  in  equation  (102)  are  quadratic  and  linear 
forms  in  arbitrary  correlated  Gaussian  RVs,  the  required  joint 
MGF  fjz(  2)  is  given  in  reference  1,  equation  (23),  by  setting 
M  *  2.  In  particular,  let  RV 

z2  =  x'  C  x  +  v'  x  +  u  ,  (109) 

where  square  matrix  C  is  positive  definite.  Then,  expressing 

z2  =  (x  +  C-1  v/2 ) '  C  (x  +  C-1  v/2 )  +  u  -  j  v'  C-1  v  ,  (110) 

it  can  be  seen  that  no  matter  what  value  random  vector  x  takes 
on,  random  variable  z2  must  satisfy 
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V  . 


(Ill) 


Z2  >  U  -  J  v'  C 

If  the  right-hand  side  of  equation  (111)  is  non-negative,  then 
the  results  in  equations  (102)  through  (108)  are  directly 
applicable,  if  rv  x  is  unable  to  take  on  the  vector  value 
-C  1  v/2 ,  then  the  non-negative  restriction  on  the  right-hand 
side  of  equation  (111)  can  be  eased.  The  general  case  of 
arbitrary  z2  values  in  ratio  (102)  is  undertaken  in  the  appendix. 
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SUMMARY 


A  method  for  obtaining  SPAs  for  some  M-dimensional 
probabilistic  problems  of  great  practical  interest  has  been 
derived  and  then  applied  to  some  representative  examples.  The 
accuracy  of  the  SPAs  has  been  verified  by  using  a  combination  of 
simulation,  asymptotic  forms,  and  Gauss-Hermite  quadrature  in  M 
dimensions.  The  SPA  with  a  first-order  correction  term  has  been 
the  model  for  these  numerical  results. 

The  main  relation  utilized  is  given  by 

f  du  p  (u)  g(u)  =  — - — rr  F  dX  //  (  X)  y(X)  ,  (112) 

J  Z  ( i2  Jt)W  J  2 

where  joint  PDF  p  (u)  and  joint  MGF  /y  (X)  are  a  Laplace  transform 

pair,  as  are  functions  g(u)  and  y(X).  M-dimensional  relation 

(112)  is  exact;  however,  numerical  evaluation  of  either  side 

requires  adoption  of  an  approximation  technique.  In  particular, 

the  approach  adopted  is  to  utilize  SPAs  with  a  correction  term; 

this,  in  turn,  requires  that  fourth-order  partial  derivatives  be 

evaluated  for  the  joint  CGF  X„(X)  =  In  // _ ( X )  . 

z  z 

For  overlapping  positive  functions  pz(u)  and  g(u),  there  is 

only  one  SP  on  the  real  axes  in  the  common  ROA  of  the  product 

fj  (X)  y(X)  in  the  transform  domain.  This  has  been  the  experience 
z 

in  the  numerical  results  here.  However,  if  p  (u)  or  g(u)  are  not 

z 

positive,  the  SPs  can  lie  anywhere  in  the  complex  X  planes. 
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All  the  g(u)  -  g(u1# . . . ,uM)  functions  considered  here  were 
separable  in  their  arguments  (u1#...fuM),  except  for  one 
two-dimensional  case.  This  separability  enables  the  use  of 
one— dimensional  Laplace  transform  tables  in  order  to  obtain  the 
corresponding  y(X)  function.  If  g(u)  is  not  separable, 
determination  of  y(X)  can  be  extremely  difficult  and  a  definite 
impediment  to  progress. 

Considerable  storage  and  execution  time  can  be  required  for 
some  of  these  M-dimensional  problems.  For  example,  the 
determination  of  the  first-order  PDF  of  the  maximum  of  a  set  of  M 
dependent  RVs  requires  that  M  M-dimensional  integrals  and  SPs  be 
evaluated.  Also,  some  combinatorial  problems  quickly  encounter 
large  increases  in  the  number  of  alternatives  that  must  be 
considered  for  an  exact  representation  of  the  probability  under 
investigation.  These  limitations  serve  to  restrict  the  use  of 
the  technique  in  practice. 
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APPENDIX  —  STATISTICS  OF  RATIO 


The  statistics  of  ratio  w  =  z^/z2,  as  in  equation  (102),  are 
of  interest,  where  //z(  a^,a2)  is  the  joint  MGF  of  RV  z  =  [  z^  z2J'. 
However,  now,  the  denominator  RV  z2  can  be  positive  or  negative. 
The  first-order  CDF  of  ratio  w  is 


cw(w)  =  Pr(z1  <  w  z2,  z2  >  0)  +  Pr^  >  w  z2,  z2  <  0) 


=  Pr(v  <  0,  z2  >  0)  +  Pr ( v  >  0,  z2  <  0)  = 


-  777-2  J  dXl  J  dX2  "vZ,(Xl'X2»/(-XlX2»  + 

(i2n)  i  r  2 

C1  ur 


+  77772  I  dXl  J  dX2  "vz  «Xl'X2>/(-XlX2»  • 

(x2n)  i  i  2 


( A-l ) 


The  RV  v  *  z1  -  w  z2,  as  in  equation  (104).  The  required  joint 
MGF  in  equation  (A-l)  is  given  by 

pvz  (Xi'X2}  =  E  exp[Xi(zi  "  w  *2)  +  X2  z2~\  =  //z(Xl,X2"wXl)  (A_2) 
2 

in  terms  of  the  joint  MGF  RV  z*  course,  it  is 

required  that  argument  (X^,X2~wX^)  e  ROA(^z).  Equation  (A-2)  can 
be  substituted  into  equation  (A-l)  to  get  first-order  CDF 


cw(w)  “  2  J  dXl  J  dX2  ^z^Xl,X2  wXl*^  X1X2  ^  + 

(  1 2  It )  p 

L1  cr 


+  (i2n)2  ^  dXl  ^  dX2  ^z(Xl'X2"wXl)/(  X1X2 ^  * 
Cr  C1 


(A- 3 ) 


A-l 


Thus,  two  two-dimensional  integrals  need  to  be  evaluated  in  order 
to  get  CDF  cw(w)  of  the  ratio  w  =  z^/^;  equation  (A-3)  is  the 
main  result  of  this  appendix.  The  A  function  of  equation  (48)  is 

A(A1,A2)  =  Xz(X1,X2-wX1)  -  ln(-X1X2)  (A-4) 

for  (X1,X2-wX1)  e  ROA (//2). 


If  RV  z2  is  always  positive,  then  joint  PDF  pz(z1,z2)  =  0  for 
argument  z2  <  0.  Then,  joint  MGF 

00  CO 

//z(Xl'X2)  =  J  dzl  exP(x1z1)  J  dz2  exp(X2z2)  pz(z1,z2)  (A-5) 

-®  0 

tends  to  zero  as  Re(X2)  ->  -®.  Then,  by  moving  the  two  X2 
contours  in  equation  (A-3)  toward  -®,  there  follows, 
respectively. 


A- 2 


As  an  application  of  this  approach  by  means  of  joint  MGF 
functions,  the  first-order  PDF  of  Student's  ratio  will  be 
derived.  Let  g  and  (g^,...,g  )  be  independent,  identically 
distributed  Gaussian  RVs  of  zero  mean  and  unit  variance  (without 
loss  of  generality).  Also,  let 


s 


5s 


(A-9) 


Student's  ratio  is 


(A— 10 ) 


RV  z2  is  obviously  positive,  thereby  allowing  use  of  equation 
(108).  The  joint  MGF  of  RV  z  =  [ z 2J'  is 

fjz(\lf\2)  =  Ez  exp(X1z1  +  X2z2)  =  Egg  exp ^X^  g  s  +  X2  s2J  - 

"  Es[eXP(X2®2)  exP(lXls2)]  "  i1  -  (2X2+Xl)/n]  ^  (A'11) 

2 

for  Re(2X2  +  X^)  <  n.  A  partial  derivative  with  respect  to  X2 
yields 

<42)<VX2>  -  t1  -  [2^l)/n]-Ha-1 .  (A-12) 

Therefore, 
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At  this  point,  make  the  change  of  variable  X  =  w  +  iu,  where 
u  is  real.  Then,  equations  (108)  and  (A-13)  yield,  for  all  w,u. 


pw(w)  "  Ili  J  dX  i1  +  (2w\-X2)/nj 


2,  , 


=  Tn  J  du  f1  +  ^w2  +  u2)/n] 


-Jjn-l 


( A-14 ) 


_ mjn+in _ 

(nil)35  r(Jsn)  [1  +  w2/n]Js(n+1) 


for  all  w  . 


(A-15 ) 


The  last  integral  on  u  was  determined  from  reference  4,  equation 
3.241  4.  Exact  result  (A-15)  is  Student's  PDF  for  ratio  w  in 
equations  (A-9)  and  (A-10). 


Alternatively,  the  SPA  to  equation  (A-14)  is  obtained  by 
observing  that 


A(X)  *  -(J$n+1)  ln[l  +  ( 2wX-X2  )/n]  , 


dA(  X) 
dX 


(n+2) (X-w) 
n+2wX-X2 


and 


d2A( X)  _  .  n+w2+( X-w) 2 

dX2  ("+2)  [n+w2-(X-w)2]2 


( A- 1 6 ) 


( A- 1 7 ) 


The  real  SP  is  obtained  from  equation  (A-16)  as  X  =  w,  for  which 
^2  *  2 

X  -  2wX  =  -w  <  n,  as  required  in  equation  (A-13).  Then, 
equation  (A-17)  yields  A2  =  (n+2)/(n+w2)  at  the  SP  X  =  w,  while 
equation  (55)  yields  SPA 
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for  all  w  . 


( A-18 ) 


This  SPA  ag  has  precisely  the  same  variation  with  argument  w  as 

does  exact  result  pw(w)  in  equation  (A-15);  however,  the  ratio 

k 

a0/pw(w)  starts  at  (n/6)^  -  0.7236  for  n  =  1  and  increases 
towards  1  as  n  increases.  The  first-order  corrected  SPA  ag  is 
ag  exp( 0 .75/( n+2 ) ) ,  which  is  0.9291  of  pw(w)  in  equation  (A-15) 
for  n  =  1. 

In  general,  the  A(X)  function  of  equation  (108)  is 

A(  X)  =  In  yu2(X,  -wX)  ,  (A-19) 

where  a  subscript  notation  for  PDs  has  been  adopted.  But,  since 
//(a,  0)  =  exp[ x(  a,  |3)  ] ,  there  follows 

^2  (  a,  0 )  =  A  (a,  0 )  X2(a,|S)  H  //(a,S)  <|»(a,|3)  ,  (A-20) 

leading  to 

A(X)  =  X(X,  -wX)  +  In  <|>(X,  -wX)  =  X  +  In  <(>  .  (A-21) 

Therefore , 

A'(X)  -  xx  -  W  x2  +  (*!  -  W  4*2 )/+  ( A-22 ) 

and 

A"  ( X )  =  xn  -  2  w  x12  +  w2  x22  + 

+  (*n  -  2  w  4>12  +  w2  ♦22)/+  -  (4X  -  w  <(>2 ) 2/<f>2  .  ( A-23 ) 
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Since  $  =  *( X,  -wX)  =  X2 (X,  -wX) ,  the  terms  in  equations  (A-22) 
and  (A-23)  are  explicitly 


*1  X12  '  *2  w  X22  ' 

and 

*11  =  X211  '  *12  =  -w  X221  '  *22  =  w  X222  * 


This  information  is  sufficient  to  locate  the  SP  of  A(X)  in 
equation  (A-19)  and  determine  SPA  aQ  to  equation  (108). 
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